
#include <iostream>
#include <fstream>
#include <complex>

int main()
{
    std::ofstream fout1("test1.dat");  // Flat Sky
    std::ofstream fout2("test2.dat");  // Spherical near equator
    std::ofstream fout3("test3.dat");  // Spherical around N pole.

    double ra0 = 0.; // RA of center for test2.

    // units here are degrees.
    for (double x = -5.; x < 5.001; x += 0.01) {
        for (double y = -5.; y < 5.001; y += 0.01) {
            double r = sqrt(x*x+y*y);
            double theta = atan2(y,x);
            if (r == 0.) continue;
            double g = 1.e-3 / r;
            if (g > 0.5) continue;

            //
            // flat sky:
            // 
            double g1 = -g * cos(2.*theta);
            double g2 = -g * sin(2.*theta);
            fout1 << x <<"  "<< y <<"  "<< g1 <<"  "<< g2 <<std::endl;

            // 
            // Spherical near equator:
            //
            
            // Use spherical triangle with A = point, B = (0,0), C = N. pole
            // a = Pi/2
            // c = 2*atan(r/2)
            // B = Pi/2 - theta
            
            // Solve the rest of the triangle with spherical trig:
            double c = 2.*atan( (r*M_PI/180.) / 2.);
            double a = M_PI/2.;
            double B = x > 0 ? M_PI/2. - theta : theta - M_PI/2.;
            if (B < 0) B += 2.*M_PI;
            if (B > 2.*M_PI) B -= 2.*M_PI;
            double cosb = cos(a)*cos(c) + sin(a)*sin(c)*cos(B);
            double b = std::abs(cosb) < 1. ? acos(cosb) : 0.;
            double cosA = (cos(a) - cos(b)*cos(c)) / (sin(b)*sin(c));
            double A = std::abs(cosA) < 1. ? acos(cosA) : 0.;
            double cosC = (cos(c) - cos(a)*cos(b)) / (sin(a)*sin(b));
            double C = std::abs(cosC) < 1. ? acos(cosC) : 0.;

            //std::cout<<"x,y = "<<x<<','<<y<<std::endl;
            //std::cout<<"a = "<<a<<std::endl;
            //std::cout<<"b = "<<b<<std::endl;
            //std::cout<<"c = "<<c<<std::endl;
            //std::cout<<"A = "<<A<<std::endl;
            //std::cout<<"B = "<<B<<std::endl;
            //std::cout<<"C = "<<C<<std::endl;

            // Compute ra,dec from these.
            // Note: increasing x is decreasing ra.  East is left on the sky!
            double ra = ra0 + x>0 ? -C : C;
            double dec = M_PI/2. - b;
            ra *= 180. / M_PI;
            dec *= 180. / M_PI;
            //std::cout<<"ra = "<<ra<<std::endl;
            //std::cout<<"dec = "<<dec<<std::endl;

            // Rotate shear relative to local west
            std::complex<double> gamma(g1,g2);
            double beta = M_PI - (A+B);
            if (x > 0) beta = -beta;
            //std::cout<<"gamma = "<<gamma<<std::endl;
            //std::cout<<"beta = "<<beta<<std::endl;
            std::complex<double> exp2ibeta(cos(2.*beta),sin(2.*beta));
            gamma *= exp2ibeta;
            //std::cout<<"gamma => "<<gamma<<std::endl;
            fout2 << ra <<"  "<< dec <<"  "<< real(gamma) <<"  "<<imag(gamma) <<std::endl;

            //
            // Spherical around N pole
            //

            dec = 90. - c * 180./M_PI;
            ra = theta * 12. / M_PI;
            fout3 << ra <<"  "<< dec <<"  "<< g <<"  "<<0. <<std::endl;
        }
    }

    return 0;
}

